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The stability of cortical function depends critically on proper regulation. Under conditions of 
migraine and stroke a breakdown of transmembrane chemical gradients can spread through cortical 
tissue. A concomitant component of this emergent spatio-temporal pattern is a depolarization of 
cells detected as slow voltage variations. The velocity of ~ 3 mm/min indicates a contribution of 
diffusion. We propose a mechanism for spreading depolarizations (SD) that rests upon a nonlocal or 
non- instantaneous feedback in a reaction-diffusion system. Depending upon the characteristic space 
and time scales of the feedback, the propagation of cortical SD can be suppressed by shifting the 
bifurcation line, which separates the parameter regime of pulse propagation from the regime where 
a local disturbance dies out. The optimisation of this feedback is elaborated for different control 
schemes and ranges of control parameters. 



During migraine and stroke neurological symp- 
toms occur representing pathological events that 
spread through the cerebral cortex. While these 
clinical observations have been known for a long 
time, only recently direct measurements were 
made. Two studies have revealed common spatio- 
temporal wave patterns, one using functional 
magnetic resonance imaging in a migraine pa- 
tient [1] and another using electrodes placed di- 
rectly on the exposed cortical surface to record 
electrical activity in a stroke patient [2]. The 
observed spatio-temporal patterns in the cortex 
constitute examples of excitable behavior that ev- 
idently emerges from pathological pathways. Spa- 
tial systems that exhibit the emergent property 
that activity breaks away from a local stimulation 
site are called excitable media [3, 4]. The capacity 
to propagate pulses is the distinguishing feature 
of excitability in spatial systems. As a mechanism 
for shifting the onset of excitability in a reaction- 
diffusion system we propose failure of nonlocal or 
non-instantaneous feedback control. 



I. INTRODUCTION 

The propagation of pathological states is a particu- 
lar aspect within the complex bidirectional relation be- 
tween migraine and stroke Here we investigate this 
aspect, in particular, how cortical tissue when modeled 
as an excitable medium becomes susceptible to spreading 
events. Our knowledge about the mechanisms of propa- 
gation is still incomplete. It is generally believed that a 
common reaction-diffusion process, called cortical spread- 
ing depression (CSD) @, 0, H, S [lo| , captures essential 
features of the observed spreading phenomena during mi- 
graine and stroke. What makes cortical tissue suscepti- 
ble to CSD has not been determined. Since the smooth 
lissencephalic cortex of animals is much more suscepti- 
ble to CSD than the convoluted cortex of human, it was 
suggested that CSD in humans occurs very close to the 
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FIG. 1: Classical measurement of electrophysiological changes 
during cortical spreading depression Upper four traces: 
logarithmic representation of changes in Na^, , Ca^"*", 
and H'^ concentration, lower two traces: extracellular poten- 
tial shift Ve, and recording from a single neuron (single unit 
activity) vs. time. (Modified from [^) 

onset of this emergent property [ll|, [12] . 

CSD is locally characterized by a nearly complete 
temporary breakdown of ion homeostasis, i.e., a stable 
stationary state (Fig. [T|) 0. Most electrophysiological 
changes follow a similar temporal course and return to 
normal after about one minute . A prominent signal 
is the slow negative potential shift Ve during CSD. It 
led to the term spreading depolarization (SD) to classify 
this and related phenomena as depolarization waves in 
the cerebral cortex [l3| . Example of related phenomena 
are SD waves that contribute to progressive deteriora- 
tion in regions surrounding an infarct core in stroke pa- 
tients. These SD waves are called periinfarct depolariza- 
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tions (PID) We use the term SD for both cortical 

spreading depression and periinfarct depolarizations in 
this paper. In this terminology spreading depolarization 
is the more general term. The respective form of SD de- 
pends mainly on differences in the energy state of cortical 
tissue. In particular, differences between CSD and PID 
concern increased blood flow compensating the increased 
energy demands, which is typical for CSD occurring in 
healthy tissue during migraine but is reduced or missing 
in PID during stroke. 

Differences in brain regions, e.g., concerning the cy- 
toarchitecture or the distribution of ion channel types, 
also modify SD. We disregard the detailed pathophys- 
iological characterization of the process, because it is 
still incompletely understood. Instead, in this study, 
SD is modeled with a standard reaction-diffusion system 
of activator-inhibitor type. Differences are reflected in 
the choice of the parameter values of the system. We 
extend the reaction-diffusion model by different kinds 
of local feedback signals. We propose that one way to 
think of such local feedback signals is in terms of control 
[l5| . In this view, the feedback represents intrinsic corti- 
cal control mechanisms that reduce cortical susceptibility 
for CSD by stabilizing the physiological state of cortical 
tissue. Consequently, their failure under certain patho- 
logical conditions leads to the emergence of CSD, e.g., 
attributed to an underlying cortical hyperexcitability in 
migraine or due to low energy levels in ischaemic 
tissue during stroke [ij]. 

In this study, different feedback mechanisms with char- 
acteristic time and space scales are considered. Possi- 
ble physiological basis of the feedback signals are mo- 
tivated and discussed. We find that the excitability of 
the reaction-diffusion system can either increase or de- 
crease. For each considered feedback scheme mainly the 
sign of the coupling strength determines whether feed- 
back can suppress wave propagation, i.e., stabilize the 
homogeneous steady state of the tissue. A specific feed- 
back scheme, in which activator and inhibitor variables 
are cross-coupled by non-local connections, indicates that 
opposed signs in the coupling strength for short-range 
and long-range connections are favorable for controlling 
the homogeneous steady state. This is a typical neuronal 
network connectivity pattern called Mexican-hat connec- 
tivity. 



II. THE MODEL 

While all evidence suggests that a reaction-diffusion 
process captures the essential features of SD, it is also 
clear that some features are oversimplified if we model 
SD with a standard reaction-diffusion mechanism of 
activator-inhibitor type. To create a more accurate al- 
beit still generic model, we extend the standard reaction- 
diffusion model by a feedback mechanism. We introduce 
the feedback signal pathway in two variations. First, we 
introduce long-range connections as a spatial coupling in 
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FIG. 2: Diagram of the spatial (r) and temporal (r) charac- 
teristics (solid line) of the feedback signal, which is used to ex- 
tend the FitzHugh-Nagumo (FHN) reaction-diffusion model. 
One of the two FHN variables is fed back via connections 
with length 5. The length is defined by the radial coordinate 
r of a local polar coordinate system, which lies parallel to the 
cortical surface and has its origin located where the connec- 
tion terminates. Due to a finite propagation velocity Cs and 
a latency time n of the feedback signal, a total time delay 
r — Ti + c~^5 occurs. A kernel function k{r) describes the 
spatial distribution of the connections, here exemplarily rep- 
resented by a Gaussian (dashed). Two physiologically moti- 
vated limit cases of this general feedback type are considered: 
▲ instant feedback along long-range connections, and ■ local 
time-delayed feedback. 

addition to diffusion. This is motivated by the fact that 
high-frequency spikes in the population activity were ob- 
served up to millimeters ahead of the approaching front 
of SD [13, Ell- Second, we consider a localized time de- 
layed feedback in the cortex, which can account for a 
slow feedback signal originating from the neurovascular 
coupling. Mathematically, these two extensions can be 
viewed as two limit cases of a general signal pathway 
(Fig. m . Such a pathway serves as an additional coupling 
mechanism that is not included in a standard reaction- 
diffusion model. With this approach, we aim to describe 
universal features of excitable media modified by a feed- 
back loop. 



A. Reaction-Diffusion IVlodel 

We model the propagation of an initially localized 
breakdown of ion homeostasis in cortical tissue (Fig. [T]) by 
an activator u and an inhibitor v as dynamic variables. 
They are coupled by kinetic reaction rates /(m, v) and 
g{u,v). We assume that only the activator species dif- 
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fuses in the medium. The equations are 

^ - f{u,v) + DV'u (1) 

-Ql = ^9M. (2) 

The parameter e is the time scale ratio in the local dy- 
namics of activator and inhibitor variables. The local 
spatial coupling is introduced by the diffusion term in 
Eq. ([T]) with the diffusion coefficient D. Without loss of 
generality, the value of D is normalized to unity. The 
choice of this value only scales the spatial coordinate. 

In principle, there are two approaches to obtain the 
reaction rates f{u,v) and g{u,v). In a bottom-up ap- 
proach, the reaction rates must be derived from a mi- 
croscopic biophysical model of SD, where all major elec- 
trophysiological properties are represented. There is no 
consent so far on the mechanism of SD. Some models 
aim to provide a complete but rather local description 
of SD [l^ [23, [2l|, that is, for processes occurring in a 
single neuron and its surrounding compartments. These 
models often introduce more than 20 dynamical vari- 
ables. Therefore, a bottom-up approach with a sub- 
sequent reduction to two major activator and inhibitor 
agents seems not to be amenable. Nevertheless, it was 
shown that a two-species activator-inhibitor model can 
reproduce the spatio-temporal pattern of SD [ll|, [l^, ■ 
Even a reaction-diffusion model with a single species was 
successfully introduced for the purpose of an order-of- 
magnitude estimate of the expected propagation velocity 
of SD 0, § . Such calculations follow essentially a top- 
down approach, which we also adopt here. The basic 
activator-inhibitor mechanism, which we describe in the 
next paragraph, can be further expanded (top-down) to 
facilitate detailed investigation of further pathways and 
variables relevant to the study of specific questions con- 
cerning SD. Successful examples of such an approach are 
computational studies [H, 0] that supported a contro- 
versial hypothesis, namely that cortical tissue surround- 
ing an infarct core dies because of the metabolic stress 
imposed by multiple SD waves. 

The variables u and v assume the roles of activator and 
inhibitor, respectively. Their kinetic functions f{u,v) 
and g{u, v) are given by a cubic nonlinearity and a linear 
function, respectively, 

f(u,v) = a(u - y) - -y (3) 

g{u,v) — u — f3 — "fv (4) 

with parameters a, (3, and 7. This is the FitzHugh- 
Nagumo (FHN) system, which is widely used as a generic 
model of excitable media [i^ll^,!!^. Before we introduce 
the feedback, we compare the simulated spatial-temporal 
patterns with the ones observed during SD. In particu- 
lar, the pulse profile and velocity in the FHN system are 
related to the corresponding quantities in SD (Fig. [1]). 
Following this approach, we can estimate realistic values 
of the parameters of the FHN system. 



B. FitzHugh-Nagumo system with feedback 

To extend the standard FHN reaction-diffusion model 
we assume that a nonlocal feedback signal s{x,y,t) is 
coupled back into the medium at any point (x, y) as 

s{x,y,t) ^ K / k{r) (5) 

JQ Jo 

(w{x + Xr, y + yr, t — t) ~ w{x, y, t)) dr d4> 

where r and (j) are the radius and the azimuthal angle, 
respectively, of a local polar coordinate system in (x, y) 
that lies parallel to the cortical surface, i.e.,Xr = r cos (j), 
and yr = rsin(f). The variable w may be chosen as ei- 
ther the activator u or the inhibitor v. The function k{r) 
is the kernel function describing the spatial distribution 
of the pathway, which is typically peaked at a distance 
r — S, and K is the coupling strength. The parameter 
T = Tp{r) + Ti, with Tp{r) = r/cs, is a time delay com- 
posed of a latency t; and a propagation delay Tp. The 
latter is due to the finite signal propagation velocity Cg. 
A schematic diagram of the signal s{x,y,t) is shown in 
Fig. [21 The feedback is, on the one hand, characterized 
by the parameter r, and the kernel function fc(r), i. e., S. 
We call this characterization the type of coupling of the 
feedback signal. On the other hand, there are two choices 
of w, namely the activator u and the inhibitor v. Further- 
more, each type can either be fed back to the activator 
u or inhibitor v rate equation, i.e., Eq. ([T]) or Eq. ([2]), 
respectively. This allows for four different combinations, 
named schemes, two self-coupling and two cross-coupling 
schemes for each type of coupling. The two self-coupling 
schemes are referred to as "mm" if w = u and "mi;" if 
w — V, and the two cross-coupling schemes are referred 
to as "mm" if w — u and "mm" if w = m. For example, in 
coupling scheme uu the activator is used to compose the 
feedback signal {w = u), which in turn is included in the 
activator rate equation (i.e., self-coupling). In a vector 
short-hand notation, the FHN system 

| = m withc^O (6) 

is replaced by 

^^F{0 + a{x,y,t), (7) 
where the coupling is given by 

with w — u or w — V. 

Note that we have constructed the nonlocal feedback 
signal ([5]) such that it contains the difference between 
the remote and local values of w. As a consequence, 
the signal s(x,y,t) tends to zero if the pattern is ho- 
mogeneous and stationary. The motivation for such a 
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FIG. 3: Scheme of the cortical surface during cortical spread- 
ing depression. Red: blood vessels formed by endothelial 
cells. The cerebral blood flow starts at arterioles and leads to 
smaller capillaries, forming the capillary bed. Blue: neurons, 
green: glial cells. A pulse of cortical spreading depression 
moves from right to left, indicated by the shadowed region. 
The extracellular concentration changes of potassium [A'"'"]e, 
the extracellular potential shift Ve, and the single unit activ- 
ity are schematically related to the pulse on the lateral front 
surface. Note that the size of the different cells types are en- 
larged by about a factor of 10 in relation to these measured 
signals. The two considered feedback mechansim are also il- 
lustrated: (i) long-range instantaneous lateral connections of 
length 5 between neurons (shown on the lateral front surface), 
(ii) a time delayed feedback control (r) of Pyragas type due 
to local response of the neurovascluar unit (indicated on the 
right lateral side surface). 



feedback is that this choice eliminates the feedback sig- 
nal in the case of successful suppression of wave prop- 
agation. The signal s{x,y,t) can be viewed as an in- 
trinsic noninvasive control because it preserves the ho- 
mogeneous steady state as the fixed point of the uncon- 
trolled {K — 0) FHN system. In this way, it can be 
compared with the common time-delayed feedback con- 
trol method introduced by Pyragas [2g|. This method 
has been widely used with great success in problems in 
physics, chemistry, biology, and medicine [l5l| including 
reaction-diffusion systems (2^, 113, HH, [H, l33l|. In particu- 
lar, it was demonstrated that it can be used to control the 
coherence and the timescales of noise-induced oscillations 
in a single FHN system |34|. ISSl i36i] and in two coupled 
excitable FHN systems [33, |3^ as well as noise-induced 
patterns in reaction-diffusion systems f39|, HO, |4l|, \^ and 
wave propagation in excitable media [43^] . This motivates 
our efforts to investigate whether a failure of such an in- 
trinsic noninvasive control scheme can explain the onset 
of excitation spread in a spatially continuous FHN sys- 
tems as a model of spreading pathological processes in 
the cortex during migraine and stroke. 



1. Instant feedback along long-range connections 

The first type of coupling that we consider is a limit 
case of Eq. ([5]). It is motivated by the observation of 
increased neuronal activity millimeters ahead of the ap- 
proaching SD front. Although the cause and the ef- 
fect of this activity remains unclear, it can be assumed 
that such activity is induced by SD via lateral cortical 
connectivity patterns (Fig. In this case the signal 
propagation velocity Cs is several orders of magnitude 
faster than the velocity c of SD, i. e., Cg ^ c (SOyitm s~^ 
in physical units). Consequently, the propagation delay 
Tp can be neglected. Furthermore, the time scale of typi- 
cal changes in electrophysiological activity is much faster 
than the time scale of changes due to the breakdown of 
ion homeostasis. The latter time scale is of the order 
of seconds (Fig. [1]), which is also in agreement with our 
simulations of the FHN system if we dimensionalize the 
model equations (see section HI below) , while the former 
time scale, i.e., the delay due to chemical or electrical 
transmission at synapses, is of the order of milliseconds. 
Therefore, it is reasonable to assume that the latency ti 
of s{x,y,t) is much smaller than the time scale of SD 
and that t; can also be neglected. This assumption can 
fail for metabotropic ion channels, like metabotropic glu- 
tamate receptors, which have increased open probabili- 
ties in the range of seconds after their activation. How- 
ever, for ionotropically mediated activity we assume that 
Tp « and T; « 0, and therefore that the additional 
coupling signal s{x, y, t) is, in this case, an instantaneous 
process (r — 0). Furthermore, we assume that the con- 
nectivity pattern is rather localized around the typical 
coupling length (5, and therefore, restricting ourselves to 
one spatial dimension, x, we approximate the kernel func- 
tion by ^-functions k{r) w 5{x — 5) + 5{x + 5). Taking 
these limits, Eq. ^ reduces to 

ss{x, t) = K(w{x -6,t)- 2w(x, t) + w{x + S, t)) (9) 

Note that Eq. ([9]) only in the limit (5^0 becomes a stan- 
dard diffusion term, however, for large S, as considered 
here, it adds a novel type of nonlocal spatial coupling to 
the FHN model. 



2. Local time-delayed feedback 

The second type is also a limit case of Eq. ([5]). We 
consider a scenario that takes into account the coupling 
within the neurovascular unit. Therefore, it is partic- 
ularly important for the dynamics of periinfarct depo- 
larizations during stroke progression. The cerebral blood 
flow (CBF) is tightly regulated by the neurovascular unit 
to meet the brain's metabolic demands. These demands 
are extremely high during SD. The dynamics of the neu- 
rovascular unit are governed by various different cells 
types (Fig. [3]) . The interaction takes place between en- 
dothelial cells building the vessel walls, mural cells con- 
trolling vessel diameter, and glia and neurons. For the 
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purpose of our simplified scheme, that is, describing uni- 
versal features of reaction-diffusion models of SD modi- 
fied by an additional signal pathway, we concentrate on 
some key features of blood flow regulation that occur 
within the time and space scale given by a single passage 
of SD and that mimic the influence of the neurovascular 
unit. 

The local CBF starts at arterioles, that is, blood vessels 
that extend and branch out from an artery and lead to 
capillaries. Capillaries are the smallest vessels, measur- 
ing 5 — 10/xm in diameter. They form the capillary bed, 
a local network supplying the brain tissue. The CBF is 
regulated by two types of mural cells. Smooth muscle 
cells regulate arterioles, and pericytes regulate capillar- 
ies. The majority of the innervation of cerebral blood ves- 
sels terminate near capillaries suggesting that blood flow 
control is initiated in capillaries [44| . Therefore, we need 
to consider that increased activity during SD produces 
an initially localized hemodynamic response evoked by 
pericytes. To describe the action of the neurovascular 
unit on SD, the temporal and spatial distribution of the 
hemodynamic response needs to be considered. 

First, we consider the spatial distribution of the hemo- 
dynamic response, which is determined by the vascular 
architecture. A precise spatial coordination of segmental 
vascular resistance is needed to effectively increase blood 
flow in a larger cortical area. Propagated vascular re- 
sponse signals are utilized to achieved this. It was shown 
that after the pericytes evoke a local capillary constric- 
tion, a pulse of constriction propagates at about 2/ims^^ 
to distant pericytes [ill. Furthermore, a dilatory signal 
mediated by release of vasoactive agents propagates from 
metabolically active cells and evokes a remote response 
in upstream precapillary arterioles (45j . Considering the 
slow velocity of 2/zms~^ we can assume that global ef- 
fects mediated by arterioles take place behind SD, i.e., 
after its passage, and therefore do not have a direct influ- 
ence on the propagation of the wave front having passed. 
We assume that the response of the capillaries is rather 
localized, so that we can approximate the kernel function 
k{x) ~ 6{x). Consequently there is only a local response. 
Still there exists a latency delay r; due to slow metabolic 
effects when a local hemodynamic response is evoked by 
pericytes. Taking these limits, Eq. ([S]) now reduces to 

Srix^t) ^ K{wix,t -r) -w{x,t)), (10) 

which is identical with the noninvasive time-delayed feed- 
back signal first introduced by Pyragas [1^ for chaos con- 
trol. This control method will now be applied to travel- 
ling pulses in excitable media. The latency delay time r 
will be varied in our simulations. 

Note that during stroke progression, the recruitment 
of periinfarct tissue into the infarct core is probably me- 
diated by wave trains of periinfarct depolarizations [l^ . 
Such periodic wave patterns are likely to be influenced 
by large scale hemodynamic responses. After the leading 
SD has initiated a slow propagating wave of constriction 
[43 | , subsequent waves in the wave train formation can be 



influenced by this. This is beyond the scope of our study. 
To model a single SD pulse, we include the neurovascular 
unit by a localized time-delayed feedback signal Sr{x,t) 
evoked from the change in blood flow which in turn is 
caused by pericytes in response to changes in neural ac- 
tivity. 



III. RESULTS 

It was suggested that SD in humans occurs close to the 
bifurcation of the onset of spreading activity in excitable 
media 0,[l3|- Therefore, we first determine the location 
of this onset in the FHN system. Then we estimate the 
parameters of the FHN system from experimental data 
on SD (Fig. [1]) and compare them with corresponding 
quantities in brain tissue. Finally, we investigate how 
the onset of excitability is influenced by the two types 
of coupling, ss{x,t) and Sr{x,t), introduced in the previ- 
ous section, in particular, how the onset is influenced by 
the newly introduced spatial and temporal scales 6 and r, 
respectively. To achieve this, we first consider the uncon- 
trolled FHN system in the regime where sustained pulse 
propagation is possible and investigate whether or not 
the additional coupling schemes suppress pulse propaga- 
tion. In this context, we can view the coupling scheme 
as a control scheme with the control goa l to stabilize 
the homogeneous steady state [i^ 13. Second, we 
choose parameter pairs (Kq^Sq) and {Kq,tq) that sup- 
press pulse propagation and determine the location of 
the onset of propagation in the FHN system with these 
coupling schemes. Furthermore, we determine the loca- 
tion of the onset of propagation in the FHN system for 
those cases where the coupling schemes cannot suppress 
pulse propagation. 



1. Boundary of pulse propagation 

The spread of an initially local activity arises in the 
system described by Eqs. if a critical parameter 

value is crossed above which the medium is susceptible for 
sustained propagating excitation patterns ^, ,49, , ,50l . l5l| . 
In a ID medium this border is obtained by finding the 
regime where stationary solutions exist in a co-moving 
frame for specific activator and inhibitor profiles. In the 
co-moving frame ^ — x + ct, the coupled partial dif- 
ferential equations for activator and inhibitor variables 
Eq. (HI)-© transform to a system of second order ordi- 
nary differential equations for U{£,) and V{£,) 

u{x,t) = C/(0 

v{x,t) = ViO (11) 

which can be further transformed to a system of cou- 
pled first order ordinary differential equations of three 
variables, Ui^, V{^), Wi^ = dU{S,)/d£, Due to 
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FIG. 4: Parameter space of the FHN system at 7 = 0.5 
and a = 1. The dimensionalized diffusion coefficient is color- 
coded as a function of the time scale ratio e and the ex- 
citability threshold parameter p. The excitable regime (color- 
shaded) and the non-excitable regime (white) are separated 
by the boundary of propagation (thick black line). Above 
this boundary traveling pulse solutions do not exist. Inset (a) 
shows a typical pulse profile corresponding to the black dot in 
the main figure. Inset (b) shows homoclinic orbits of the FHN 
system in the co-moving frame with velocity c as a function of 
13 for different values of e. In this diagram, the location of the 
boundary of propagation is determined by the value of /3 at 
the turning point of each curve and the corresponding e value, 
exemplarily shown by cyan and magenta dots for e = 0.001 
and e = 0.2, respectively, in inset (b) and the main figure. 



this transformation, c, the propagation velocity, is intro- 
duced as an additional parameter. In this system, homo- 
clinic orbits can be found using a continuation method 
[5^. They correspond to travelling pulse solutions in 
Eq. (Il])-([2]) (Fig.m inset b). The transition into a region 
where homoclinic orbits exist marks a bifurcation of codi- 
mension one. In the parameter space of excitable media 
this bifurcation separates the regime supporting travel- 
ling pulses from the non-excitable regime. This border 
is shown in Fig. 2] in the parameter plane (e, (3) at the 
section a = 1 and 7 = 0.5. As expected, the transi- 
tion into the non-excitable regime is achieved by increas- 
ing either [3 or e. Increasing /3 increases the threshold 
of the excitable medium, while increasing e changes the 
time scale ratio between the dynamics of activator and 
inhibitor, with larger values of e corresponding to faster 
inhibitor dynamics. 



2. Estimate of space and time scales 

The main reason to introduce spatial and temporal 
units in our FHN system is to obtain a better under- 
standing of the typical values of the spatial and temporal 
coupling scales, 5 and t, respectively. Furthermore, the 
value of the diffusion coefficient D can be determined 



if we introduce dimensional units in the FHN system. 
Note that D is not a bifurcation parameter. In the non- 
dimensionalized FHN system Eqs. D can assume 

any value without changing the dynamics in a qualita- 
tive way. For example, changing D has no effect on the 
location of the propagation boundary shown in Fig. |31 
Varying D is equivalent to scaling space units by _D^2. 
Therefore, the physical value of the diffusion coefficient D 
can be determined if we match the patterns obtained in 
the FHN system with the ones observed experimentally 
during SD. 

First, we define characteristic space and time units, 
xq and respectively. This is done with respect to 
the spatio-temporal patterns of the uncontrolled {K = 
0) FHN system. These characteristic units are cho- 
sen such that the patterns calculated from the non- 
dimensionalized Eqs. match the measured ones 

of Fig. [T] Introducing dimensional space (X) and time 
(T) variables by a; = AT/xo and t = T/to, the FHN sys- 
tem Eqs. assumes its dimensional form with the 
dimensionalized diffusion coefficient D = Xq/Iq. 

Now we compare the simulated pulse width Ax and 
duration At with the measured pulse width AX and 
measured duration AT (Fig. [H AT « 20s). Using the 
typical measured velocity of SD C « 50/ims~^, we ob- 
tain AX = CAT « 0.1cm. The simulated pulse width 
Aa; and duration At are related by Ax = cAt with the 
non-dimensionalized propagation velocity c, which ap- 
pears as a parameter in Eq. pT|) (e.g. for e = 0.1 and 
/3 = 0.85 : Ax « 8.7 in Fig. H (inset a) and c w 0.81 
in Fig. 0] (inset b), hence At w 10.70). Hence we infer 
the space and time units xq = AX/ Ax ~ 115/im and 
to — AT/ At « 1.9s, respectively, which yields a diffu- 
sion coefficient D = Xq/^q ~ 70 • 10~^cm^s~^, which is 
a reasonable value. Note that the pulse width Ax, the 
velocity c and hence xq, io and D depend upon the pa- 
rameters (3 and e. 

In Fig. [4] the obtained values for the diffusion coeffi- 
cient D for each pair (e,/3) are shown in color coding. It 
ranges from 5 • 10~® — 100 • 10~^cm^s~^ in the main part 
of the parameter space. As expected, low values are ob- 
tained in regions far away from the boundary of propaga- 
tion. There the excitability of the system is high. On the 
other side, close to the propagation boundary, the values 
of D are higher. In this regime, the medium is weakly ex- 
citable or subexcitable. The correlation between the value 
of the obtained diffusion coefficients D and the excitabil- 
ity regimes requires some further clarification, which will 
be addressed in the discussion. However, note that the 
values of D in the parameter plane (e, (3) at the section 
a — I and 7 — 0.5 are in the range of expected values for 
a diffusing substance participating in the mechanism of 
SD, as discussed in the literature 0, d, [H, [M HI ■ 
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3. Control of spreading depressions by ss{x,t) and Sr(x,i) 

How does the additional coupling s{x,t) change the 
FHN system? To answer this, we perform simulations 
with a variety of FHN systems in the parameter plane 
(e,/3) at a = 1 and 7 — 0.5. The systems are chosen 
close to the propagation boundary (Fig. 3]) . We include 
both types of coupling, the instant feedback along long- 
range connections ss{x,t), described in section H.l, and 
the local time-delayed feedback St{x, t), described in sec- 
tion H.2. Both types are considered separately. For 
each type of coupling there are four coupling schemes, 
two self-coupling schemes (mm, vv) and two cross-coupling 
schemes (mu, uv). 

To obtain the influence of the coupling schemes on the 
FHN system, we start each individual simulation with a 
stable pulse profile of the reaction-diffusion system ob- 
tained without the couphng signal {K = 0). Then, for 
K ^ 0, we determined whether or not the pulse propa- 
gation is terminated. If the propagation is suppressed, 
we also determine how long the pulse can still propagate 
before it disappears. This distance defines the volume of 
tissue at risk (TAR), referring to the risk of cortical tissue 
surrounding a local pathological core of being recruited 
into the disturbed state . This TAR value is taken as 
a measure of the efficiency of the coupling scheme as a 
control method. 

First, we consider the four non-local schemes in 
ss{x,t), i.e., instant coupling along long-range connec- 
tions (Sec. H.l). Regions in which pulse propagation is 
suppressed are shown in Fig. [5] a-d. The TAR value is 
given by a color code in this regions. For the two cases 
where we introduce self-coupling (Fig. [5] a (w — u) and 
d {w = v)), mainly the sign of the coupling constant 
K determines the success of the coupling scheme. Pulse 
propagation cannot be suppressed with negative values 
of K. This can be intuitively understood, if we consider 
only the effect of the signal ss{x,t) on the homogeneous 
steady state. A small disturbance from the homogeneous 
state is destabilized by ss{x,t) UK < and stabilized 
UK > 0. Optimal values of S are obtained for w = u 
at (5 w Aa; and for w — v a,t 5 k, 0.5Ax. The picture 
changes for the two schemes with cross-coupling (Fig. [S] 
b (w = v) and c (ui = m)). Successful suppression of pulse 
propagation in these cases depends on both the coupling 
strength K and the coupling length 5. The change in 
the sign of K occurs at (5 « Aa;. li 5 < Ax, K must 
be positive (negative) if the coupling term ss(x,t) is fed 
into the inhibitor (activator) balance equation Eq. Q or 
(fT|), respectively, and vice versa for J Ax. In other words, 
the short and long range control domains have opposite 
signs of K. Optimal values of 5 in the short range control 
domain are obtained at 5 ~ 0.5Ax. 

Next, we consider the local time-delayed feedback cou- 
pling (Sec. II. 2). In each of the four coupling schemes 
control is only achieved if K is either positive or nega- 
tive (Fig. [5] e-h). For positive K, the pulse propagation 
is suppressed if the signal is coupled into the activator u. 



ss(x,t) 




-1 -0.5 0.5 1 -1 -0.5 0.5 1 



K K 

FIG. 5: Control planes for the four coupling schemes of the 
two types of feedback ss{x,t) (a-d) and ST{x,i) (e-h). (a,e): 
self-coupling of the activator signal (uu); (b,f): self-coupling of 
the inhibitor signal ivv); (c,g): cross-coupling with the feed- 
back signal composed of the activator {w = u) and fed back 
to the inhibitor rate equation (uv); (d,h): cross-coupling with 
the inverse configuration as compared to (c,g), i. e., (vu). The 
parameter values of the FHN system are e — 0.1, /3 = 0.85, 
a = 1, and 7 — 0.5. This system is close to the propagation 
boundary (see black dot in Fig. |4} . Regions in which pulse 
propagation is suppressed by the feedback signal are shown 
color coded by the tissue at risk (TAR) value, that is, the 
volume of tissue recruited into the pathological state before 
the pulse dies out. Low TAR values indicate that the feed- 
back signal is more efficient. Optimal values, i.e., low TAR 
values, for 5 and r are mainly located at about 0.5 times the 
activator pulse width and duration, respectively (see text). 
Red and green dots mark control parameter values {K, 5) and 
{K, t) for which the shift in the pulse propagation boundary 
is calculated in Fig. |6l 



The two domains in the (K, r) plane, where pulse prop- 
agation is suppressed, are of comparable size (Fig. O e 
(self-coupling) and f (cross-coupling)). For negative K, 
the pulse propagation is suppressed if the signal is cou- 
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pled into the inhibitor v. In this case, the cross-couphng 
scheme (Fig. [5] g) has a domain of comparable size and 
shape as for cases shown in Fig.[5]e,f. The domain for the 
self-coupling scheme is much smaller with only large val- 
ues of TAR, i. e., less efficient suppression of pulse prop- 
agation (Fig. Oh). For all four cases, the optimal delay 
time is r « 0.5Ai. 



4. Shift in tlie onset of excitability by ss{x,t) and ST{x,t) 

We now investigate the shift in the location of the prop- 
agation boundary (Fig. |4]) caused by the four coupling 
schemes of both types of coupling ss{x,t) and ST-(x,t). 
This boundary separates the excitable regime from the 
non-excitable regime. Intuitively it is clear that the sup- 
pression of pulses shifts the propagation boundary to- 
wards the excitable regime of the system with K = 0. 
If we consider the coupling schemes as control mecha- 
nisms, such a shift manifests itself as a successful control 
strategy that reduces the excitability of the system. 

First, we consider the non-local type of coupling 
ss{x,t). The parameter values for ss{x,t) are set to 
K = 0.2 and K = -0.2, and S = 0.5Ax. For each of 
the four coupling schemes with 5 = 0.5Ax only one value 
of K, either K = 0.2 or K = —0.2, suppresses pulse 
propagation (marked by green dots in Fig.E]). With this 
choice and depending upon the four coupling schemes, 
the control signal ss{x,t) lies either close to optimal 6 
values or on the opposite side, mirrored along the K = 
axis (marked by red dots in Fig. [5]), a location where 
the control goal is not achieved. In the cases of the two 
cross-coupling schemes, the value K = —0.2 for w — v 
(Fig. [5]b) and the value K = 0.2 for w = u (Fig. [5]c) lies 
in the center of the control domain with the short range 
connections. If we consider the four cases in which pulses 
are suppressed, the location of the propagation bound- 
ary is shifted into the excitable regime, indicated by dif- 
ferently colored regions in Fig. [S]a). It is evident that 
the uv and vu cross coupling schemes are most efficient. 
Note that for the least efficient coupling scheme (uu) the 
propagation boundary (edge of the green region in Fig. [5] 
a) is actually not shifted beyond the point e = 0.1 and 
/3 — 0.85 (black dot). This seems to be contradictory 
with the control domain shown in Fig. [5] a) but is consis- 
tent with the definition of the propagation boundary as 
a bifurcation line. Beyond this bifurcation line there do 
not exist any traveling pulse solutions, whereas Fig. [5] a) 
only states that a specific traveling pulse solution with a 
specific amplitude (the one of the uncontrolled system) is 
suppressed. This example illustrates the different nature 
of the information shown in Fig. [5] and Fig. [SI In the 
four cases with the opposite sign of the parameter value 
K, i. e., where pulses are not suppressed, the propagation 
boundary is shifted towards the non-excitable regime of 
the uncontrolled {K = 0) FHN system (dashed lines in 
Fig.Ela). 

Second, we consider the local, time-delayed type of 




e 

FIG. 6: Parameter space (e, /3) as in Fig. [4] showing excitable 
and non-excitable regimes for (a) non-local, (b) local time- 
delayed coupling. For the uncontrolled system {K = 0) the 
border between the excitable regime and the non-excitable 
regimes is between the black and white region. The other 
colored regions (green, red, blue) visualize the reduction of 
the excitable regime by the different coupling schemes as in- 
dicated by the location of the green dots in Fig. [5}, i.e. K 
is chosen as K — 0.2 (uu, uv, vv in (a); uu, vu in (b)), 
or K — —0.2 (vu in (a); uv, vv in (b)). For each coupling 
scheme with opposite sign of K (red dots in Fig. [SJ the ex- 
citable regime increases as indicated by the dashed lines. 



coupling Sr{x,t). The parameter values for ST-{x,t) are 
set to K = 0.2 and K = -0.2, and r = 0.5At. To 
provide a comparison with the former type of coupling 
ss{x,t), we also choose the short range domain (as for 
ss{x,t) in Fig. [5]b,c) as the reference point and use these 
values of K also for St{x, t). As before they are marked 
by green dots in Fig. [5] e,f,g if control is successful and 
red otherwise. Only for the scheme (vv) (Fig. [5] h) both 
types of coupling for K — —0.2 and K — 0.2 lie outside 
the control domain. Since the type with K = —0.2 is 
closer to a control domain, this is marked green, and cor- 
respondingly its propagation boundary is slightly pulled 
back as indicated by the red region in Fig. [6] b. Again 
we find a similar pattern as for the non-local type of 
coupling ss{x,t). The coupling schemes that suppress 
pulse propagation reduce excitability (colored regions in 
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Fig. [6] b) and the other ones increase excitability (dashed 
hnes). Note that also for the coupling scheme (vv) this 
pattern occurs, with the only difference that the prop- 
agation boundary (edge of the red region in Fig. [5] b) 
is not shifted beyond the point e = 0.1 and (3 — 0.85 
(black dot), which is in agreement with the control do- 
main shown in Fig. [5] g. 



DISCUSSION AND CONCLUSIONS 

We have shown that a weakly excitable reaction- 
diffusion system can be shifted across the boundary of 
pulse propagation (Fig. [5]) by various types of feedback 
and coupling schemes. In particular, we use non-local 
feedback ssix,t) and time-delayed feedback Srixjt). On 
the one hand, this is motivated by physiological consid- 
erations (Fig. [3]) . On the other hand, these two distin- 
guished types allow for a clear separation of the effect of 
the introduced space scale S and time scale t upon the 
spatio-temporal pattern of SD. We have shown that these 
scales 6 and r are closely related. The reason is, of cource, 
that the feedback is applied to control of travelling pulses. 
The velocity of SD relates space to time scales. Their nu- 
merical values can be compared after normalizing them 
with the pulse width Ax and pulse duration At, respec- 
tively. The common effect of both types of feedback on 
the traveling pulse is then reflected in the similar location 
of optimal control around half the pulse width and pulse 
duration for ss{x,t) and Sr{x,t), respectively (Fig. [5]). 
Therefore, the feedback mechanism, as given in a gen- 
eral form in Eq. (O, and illustrated in Fig.[2l provides a 
generic mechanism to modulate excitability in reaction- 
diffusion systems even if we consider only its two limit 
cases. 

The reaction-diffusion system we have used to demon- 
strate the modulation of excitability by feedback is the 
FHN model. It represents excitable media exhibiting ex- 
citability of type II. Furthermore, we have constrained 
our investigations to the section a = I and 7 = 0.5 in the 
four-dimensional parameter space of this model. Both 
restrictions, excitability of type II and this choice within 
the parameter space, will now be discussed. Type II ex- 
citability is related to a Hopf bifurcation, which is char- 
acterized by the fact that a periodic state bifurcates with 
a nonzero frequency [s^ . Within an exponentially small 
range of either or e after the Hopf bifurcation a tran- 
sition occurs, called canard explosion, from a state with 
small amplitude oscillations to a large amplitude relax- 
ation oscillation. The canard explosion produces the 
threshold behavior (all-or-none) of the system needed 
to exhibit excitable behavior. Another typical bifurca- 
tion scenario leading to excitabilit y is that the periodic 
state appears at a zero frequency [STI fssj , called type I 
excitability. This can be modeled by a saddle-node bi- 
furcation on a limit cycle or saddle-node infinite period 
bifurcation (SNIPER) 

It is not clear which type of excitability is the ba- 



sis of SD. However, experiments with retinal spreading 
depression indicate excitability of type II behavior [60j . 
The extracellular potassium concentration ([X+Je) was 
increased stepwise. Above a usually well preserved ceil- 
ing level of [ii^^]e = lOmM 6l| a rather spontaneous on- 
set of oscillations with finite period was observed. Due 
to the presence of noise and a rather large step width of 
A[iir+]e = 2mM, a definite conclusion is not yet possible. 
However, we assume that our results hold in principle also 
for excitable media in which each local element is of ex- 
citability type I. The difference is manifest only close to 
the bifurcating state of the individual excitable element. 
In an excitable media, that is in a spatially extended sys- 
tem, excitability is defined by the distinguishing feature 
that the excited state breaks away from a local stimu- 
lated area due to transport, which usually involves dif- 
fusion. The bifurcating state of the individual elements 
seems therefore less important than the saddle-node bi- 
furcation that leads to the emergence of traveling pulses 
(Fig.H inset b). 

The FHN system is the most generic type that shows 
excitability of type II, because the activator equation, 
Eq. ([T]), has a cubic nonlinearity, which is generic for 
bistability [g^, [g^ . Since we aim to describe generic fea- 
tures, let us briefly comment on the choice of the two 
parameters that we have fixed (a = 1 and 7 = 0.5). 
With this choice, the value of the diffusion coefficient in 
the main part of the remaining section of the parameter 
space is between 5 • 10~^cm^s~^ and 100 • 10~^cm^s~^ 
(Fig. [4|). The diffusion coefficient of [-ft^'^'Je, a substance 
often related to the activator or at least to the species by 
which SD propagates, is about 20 • 10~^cm^s~^ in aque- 
ous solution [64|. There are two counteracting effects 
that can significantly change this value in cortical tissue. 
First, due to the porous geometrical structure and small 
volume fraction of the exctracellular space, which is sim- 
ilar to a soap phase, the apparent diffusion coefficient of 
K~^]e in brain tissue is estimated to 7.2 • 10~^cm^s~^ 
11]. Second, the possibility that [K~^]e enters glial cells 
in one place and leaves elsewhere provides a form of fa- 
cilitated diffusion that is estimated to be up to five times 
more important than diffusion in the extracellular space 
[8| . In summary, a generic model that simulates traveling 
pulses with a diffusion coefficient D of reasonable order 
of magnitude seems to be appropriate. 

We propose that a failure of feedback provides a com- 
mon mechanism of the emergence of spreading depolar- 
izations in migraine and stroke. A recent study taking 
a complementary bottom-up approach, by describing SD 
in a detailed biophysical model, has concluded that the 
key to normal stability of cortical tissue is the effective 
regulation of [i^+]e by the neuron's Na-K ion pump and 
the glia-endothelial system This complements our 

conclusions. In our generic approach, our concern is not 
to identify the precise species behind activator, inhibitor, 
and feedback signals, but rather to describe spatial prop- 
erties of their interaction and the emergence of travelling 
pulses. On the contrary, a detailed biophysical model of 
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properties of a single neuron and its surrounding extra- 
cellular compartments can only infer statements about 
the local excitable state. Whether a local breakdown of 
ion and transmitter homeostasis recruits further tissue 
into the excitable state or remains restricted to the ini- 
tial focus is the clinically important question. This can 
only be answered in an extended system with sufficient 
spatial resolution, where a pulse propagation boundary 
can be defined (Fig.|4]). The existence of the propagation 
boundary has consequences for therapeutic approaches 
that aim at an effective reduction of excitability to hold 



the tissue below the bifurcating state and with that limit 
the tissue at risk [Sol. 
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